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ABSTRACT 

We describe a sampling method to estimate the polarized CMB signal from observed maps of the 
sky. We use a Metropolis- within- Gibbs algorithm to estimate the polarized CMB map, containing Q 
and U Stokes parameters at each pixel, and its covariance matrix. These can be used as inputs for 
cosmological analyses. The polarized sky signal is parameterized as the sum of three components: 
CMB, synchrotron emission, and thermal dust emission. The polarized Galactic components are 
modeled with spatially varying power law spectral indices for the synchrotron, and a fixed power 
law for the dust, and their component maps are estimated as by-products. We apply the method to 
simulated low resolution maps with pixels of side 7.2 degrees, using diagonal and full noise realizations 
drawn from the WMAP noise matrices. The CMB maps are recovered with goodness of fit consistent 
with errors. Computing the likelihood of the E-mode power in the maps as a function of optical 
depth to reionization, r, for fixed temperature anisotropy power, we recover r = 0.091 ± 0.019 for a 
simulation with input r = 0.1, and mean r = 0.098 averaged over 10 simulations. A 'null' simulation 
with no polarized CMB signal has maximum likelihood consistent with r = 0. The method is applied 
to the five-year WMAP data, using the K, Ka, Q and V channels. We find r = 0.090 ± 0.019, 
compared to r = 0.086 ±0.016 from the template-cleaned maps used in the primary WMAP analysis. 
The synchrotron spectral index, averaged over high signal-to-noise pixels with standard deviation 
(j{/3) < 0.25, but excluding ~ 6% of the sky masked in the Galactic plane, is —3.03 ± 0.04. This 
estimate does not vary significantly with Galactic latitude, although includes an informative prior. 
Subject headings: cosmic microwave background, cosmology: observations, methods: statistical, po- 
larization, radio continuum: ISM 



1. INTRODUCTION 



The Wilkinson Microwave Anisotropy Probe (WMAP) 
has mapped the sky in five frequency bands between 
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23 and 94 GHz. Measurements of the temperature 
anisotropy in the CMB have led to the establishment 
of the ACDM cosmological model. Anisotropics in the 
CMB polarization at large scales inform us about the 
ionization history of the universe, allow us to probe 
a possible signal from gravitational waves seeded early 
in the universe, and lead to improved constraints on 
cosmological parameters when combined with tempera- 
ture measuremen ts. The three-year WMAP observations 
(jPage et al.ll2007l ) showed that polarized diffuse emission 
from our Galaxy dominates the primordial signal over 
much of the sky, making accurate estimation of the CMB 
signal at large angular scales challenging. 

In this paper we describe a Bayesian framework for es- 
timating the low resolution polarized CMB maps, and 
errors marginalized over possible Galactic emission. The 
goal of this approach is to determine not only the 'best' 
estimate of the microwave background polarization fiuc- 
tuations but to determine the uncertainties associated 
with foreground re moval. A simi l ar technique has also 
been developed by lEriksen et all (|2006l . I20Q7D for esti- 
mating intensity maps, and has been applied to the three- 
year WMAP temperature maps. This work complements 
the primary analysi s of the five- year WMAP polarization 
maps described in iGold et al.l (2008), which uses tem- 
plate cleaning to estimate the CMB polarization maps. 
The basic five year WMAP results are summarized in 
iHinshaw et all (|2008[ ). This paper is structured as fol- 
lows. In 32] we describe the sampling method used to 
estimate polarized maps. In §3] we apply it to simulated 
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maps, and in §3]to the WMAP data. We conclude in ^ 

2. ESTIMATION OF POLARIZATION MAPS 

The large-scale polarized radiation observed by 
WMAP is the sum of at least three components: the pri- 
mordial CMB, synchrotron emission, and thermal dust 
emission. Here we briefly review the Galactic emission 
mechanisms, for more details see e.g.. Page et al. (2007). 
Both synchrotron and thermal dust emission are polar- 
ized due to the Galactic magnetic field, measured to 
have a coherent spiral structure parallel to the Galac- 
tic plane, as well as a significant turb ulent comp onent 
(Spitzer 1998; Beck 2001; Vallee 2005: [Haj]|[2QQ6f ). The 
effective strength of the field is of order ^10 /iG, thought 
to be split roughly e qually between the co herent and tur- 
bulent components (jCrutcher et al.l [20031 ). Synchrotron 
emission is produced by relativistic cosmic-ray electrons 
accelerated in this magnetic field (see Strong et al. (2007) 
for a review of cosmic ray propagation). For electrons 
with a power law distribution of energies 

N{E) oc E-P, (1) 

the frequency dependence of the emission is character- 
ized by antenna temperature T{h') oc with spec- 
tral index (3 = — (p + 3) /2, with typically P ^ —3 
(|Rvbicki fc Lightmanll 19791 ). However, since synchrotron 
loss is proportional to E^, older sources of electrons 
should have a lower energy distribution and a steeper 
spectral index of synchrotron emission, compared to re- 
gions of recently injected electrons. This l eads to a syn- 
chrotron index that varie s over the sky (jLawson et al.l 
119871 : I Reich fc ReichI 119881) an d is expected to s teepen 
away from the Galactic plane ([Strong et al.ll20Q7l ). with 
evidence of this behavior seen in the WMAP data 
(jBennett et al. 2003). Since the cosmic-ray electrons 
emit radiation almost perpendicular to the Galactic mag- 
netic field in which they orbit, the y can produce polar- 
ization fractions as high as ~ 75% ([Rvbicki fc LightmanI 
IT979), although integration of multiple field directions 
along a line of sight reduces this level. The fractional 
polarization observed at radio frequencies in the range 
408 MHz - 2.4 GHz is further l owered due to Fara- 
day rotation (Duncan et al.l 1 19951 : lUvaniker et al.l 1 19991 : 
IWolleben et~al...2006i ). In the WMAP 23 GHz data, the 
polarization frac tion is as hi^h as 50 % on significant por- 
tions of the sky fKogut et al]l2007[ ). 

Thermal dust intensity has been well measured 
by the IRAS and COBE m issions and extra polated 
to microwave frequencies by iFinkbeiner et al.l (|l999). 
Polarization arises since grains tend to align their 
long axes perpendicular to the Galactic magnetic 
field v ia, for example, the Da vis- Greenstein mecha- 
nism ([Davis fc GreensteinI 119511 ). and depending on 
their composition can be polarized up to a mod- 
eled maximum of ^ 20% p a rallel t o the long axes 
(e^ lHildebrand fc DragovanI ([1995r ): [D raine fc Fraissd 
([2008[ )). Observations of starlight, polarized perpendic- 
ular to the d ust grai ns, are co nsistent with this picture 
([HeilesI [200(1 BerdvugijLet ^. 2001), as are the three- 
year WMAP observations (Pa ge et al.[[2007[ : [Kogut et al.[ 
[20071 ). A population of smaller spinning dust grains 
formed of polycyclic aromatic hydrocarbons may also 
emit a significant amount of microwave radiation due 
to electric dipole rotational emission (Draine fc Lazarian 



[1999[ : [Drain e fc Lil l20(^. Thi s question is discussed 
i n e^ . flinshaw et al] ([2007[ ): [Dobler fc Finkbeined 
([2007[ ): [Gold et al.< (20Q8[ ) with respect to the inten- 
sity signal observed by WMAP. However, these small 
spinning dust gr ains are not expected to be signifi- 
cantly polarized ([Draind [2003[ ). Other mechanisms for 
producing polarized emis sion, including magnetic dust 
([Draine fc Lazarian[[l999[ ). have not been observed to be 
dominant. 

Given these polarized Galactic components, the stan- 
dard method used to clean the WMAP polarization maps 
involves subtracting synchrotron and dust template maps 
from the total, l eaving a cleaned CMB map at th e Ka, 
Q, and V bands ([Page et al.[[2007[ : [Gold et al.[[2008[ ). The 
spectral indices of the templates are not allowed to vary 
spatially, which is a sufficient approximation given the 
sensitivity of the observations. Errors are propagated 
by inflating the noise matrices to account for the un- 
certainties in the fitted coefficients of template maps 
([Page et al.[[2007[ ). In this alternative method we param- 
eterize the emission model, and use a sampling method 
to estimate the marginal mean posterior CMB Q and U 
maps in low resolution pixels. 

2.1. Bayesian estimation of sky maps 

The data, d, consist of the Q and U polarization maps 
observed by WMAP at Nc frequency channels, and is a 
vector of length 2Np x Nc. In this analysis we will use 
HEALPix A^side=8 with Np = 768 The joint posterior 
distribution for a model m can be written as 

p(m|d) oc p(d|m)p(m), (2) 

with prior distribution p(m) and Gaussian likelihood 

-21np(d|m) =x'+c, (3) 

with 

= (d - m)^N-^(d - m) (4) 

and a normalization term c. Since the noise N is uncor- 
related between channels, the likelihood can be written 
as the sum over frequency 

x' = ^[d,-m,]^N;Hd.-m,], (5) 

V 

where N^^ is the noise covariance at each channel. In this 
analysis we use WMAP low resolution maps with inverse 
noise matrices that describe the noise outsi de a process- 
ing ma sk covering ~ 6% of the sky (see e.g., Jarosik et al.[ 
([2007[ )). These masked pixels should be neglected in the 
likelihood evaluation, but to simplify numerical imple- 
mentation we include them, but set their inverse vari- 
ance to be small (less than 0.1% of the unmasked pixels' 
inverse variance). 

We parameterize the model in antenna temperature as 
the sum of three components: 

= ^afe,i.Afc, (6) 

k 

with CMB (/c=l), synchrotron emission (/c=2), and ther- 
mal dust emission (/c=3). In this analysis we will ignore 

16 The number of pixels is = 12N ^.^^, with N^i^e = 2^^^ 
where res is the 'resolution' (jGorski et al.|[2005i ). 



WMAP 5-year Polarized Map Estimation 



3 



possible polarized contributions from other components 
including spinning dust, and free- free emission. Free- free 
emission may become slightly polarized due to Thomson 
scattering by electrons in HII regions. 

The components each have amplitude vectors A/^ of 
length 2Np and diagonal coefficient matrices cxk^u of side 
2Np at each frequency. The CMB radiation is black- 
body, and we further assume that the Galactic compo- 
nents can be described with a spectral index that does 
not vary with frequency in the WMAP range. The coef- 
ficients are therefore given by 

cxi,u = f{iy% (7) 

a2,. = diag[(^/^K)^2], (8) 

cx3,u = did.g[{u/uw)^^]- (9) 

Here we have introduced two spectral index vectors /3j^ 
each of length 2A/p, pivoted at 23 GHz {i^k) and 94 GHz 
(i^w) respectively. The function /{u) converts thermo- 
dynamic to antenna temperature. We then make three 
further simplifying assumptions. First, that the spectral 
indices in Q and U are the same in a given pixel, which 
equates to assuming that the polarization angle does not 
change with frequency. Second, the spectrum of thermal 
dust is assumed to be fixed o ver the sky, with fiducial 
value = 1.7, motivated bv iFinkbeiner et all (p^999> ). 
Third, we define Ni < Np pixels within which P2 takes a 
common value, rather than allow it to take a unique value 
at each A/side=8 pixel. This is motivated by our under- 
standing of the emission process: even though we expect 
spatial variation due to the different ages of the elec- 
tron populations, the electron diffusion ra te limits how 
much the in dex can vary over a short range ([Strong et al.l 
l2QQ0l . l2QQ7r ). In this case we use iVside=2 HEALPix pixels 
(A^i=48). 

Our model m is now described by 6Np amplitude 
parameters A = (Ai,A2,A3)-^ and Ni spectral index 
parameters /3. Our main objective is to estimate the 
marginalized distribution for the CMB amplitude vector, 

p(Ai|d) = Jp{A,(3\d)dA2dAsd^, (10) 

from which we can estimate a map and covariance ma- 
trix. 

2.2. Sampling the distribution 

We cannot sample the joint distribution p(A,/3|d) di- 
rectly, so we use Markov Chain Monte Carlo methods to 
draw samples from it. It can be sliced into two condi- 
tional distributions p(A|/3,d) and p(/3|A,d), so we use 
Gibbs sampling to draw alternately from each conditional 
distribution, constructing a Markov chain with the de- 
sired joint distribution as its stationary distribution. 

We briefly review Gibbs sampling for the case of one 
A parameter and one /3 parameter: starting from an ar- 
bitrary point {Ai^Pi) in the parameter space, we draw 

(A,+i,A+l),(A,+2,A+2)... (11) 

by first drawing A^+i from p{A\/3i^d) and then drawing 
/3i-^i from p{/3\Ai-^i^ d). Then we iterate many times. 
The result is a Markov chain whose stationary distribu- 
tion is p{A, (3\d). A description of Gibbs sampling can be 



fou nd inlGelfand SmitH ([T99Qh . I Wandelt et all ([2QQ1 . 
and lEriksen et all (|2QQ7D . 

For the multivariate case A and f3 are now vectors, 
and so each vector is drawn in turn until convergence, 
producing a chain whose stationary distribution is the 
joint posterior distribution. Two distinct methods are 
used to draw the samples from each conditional distribu- 
tion, depending on whether the amplitude vector A, or 
the index vector /3, is held fixed. 

2.2.1. Sampling the amplitude vector 

For fixed /3, the conditional distribution p(A|/3,d) is 
a 67Vp-dimensional Gaussian, so one can draw a sample 
of all 6Np amplitude parameters simultaneously. The 
conditional distribution is 

p{A\(3,d)^p{d\p,A)p{A). (12) 

For a uniform prior on A, the mean. A, is found by 
minimizing 

= ^[d, - ^ ak,.AkfN-^[d, - cck,uAk] (13) 

1^ fc k 

with respect to A. This gives A = F~^x, which can be 
written in block-matrix form, 

Ai \ /Fn Fi2 Fi3 \ /xi \ 
Aa = F21 F22 F23 X2 (14) 

VF31 F32 F33/ Vx3/ 

with elements 

Fkk'=Y.al^^-^ak',u, (15) 
Xfc = ^aJ,,N;id,. (16) 

V 

Note that F is a 2Np x 2Np matrix and x is a vector 
of length 2Np. The covariance of the conditional distri- 
bution is given by F~^. In the case of diagonal noise, 
the mean and variance are estimated pixel by pixel using 
the same method. Given the mean and Fisher matrix 
of the conditional distribution, we draw a Gaussian sam- 
ple using the lower Cholesky decomposition of the Fisher 
matrix, F = LL^, with sample A^+i = A + L~^G. The 
vector G contains 2Np zero mean unit variance Gaussian 
random samples. 

For a diagonal Gaussian prior on A^^h of a/c ± cr/c, the 
expressions are modified to 

F/cfe' =F/efe' + (5fefe'Cr^^I, (17) 
X/e=X/e +Cr^^Iafc, (18) 

with posterior mean A = F~^x and variance F~^. In 
this analysis we place uniform priors on the CMB and 
synchrotron Q and U amplitudes at each pixel, but im- 
pose a Gaussian prior on the dust Stokes vector A2 = 
(Q2, U2)^ of [Q2, U2] = ± 0. 2 \fi, usin^ the dust ma p 
Irf at 94 GHz from model 8 of IFinkbeiner et all (|1999[ ). 
hereafter FDS, as a tracer of the intensity. The width of 
the prior, corre sponding to a pol arization fraction 20%, 
is motivated bv lDraine fc Fraissd ([2008), who predict the 
maximum polarization of dust grains to be about 15%. 
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Drawing this new amplitude vector is computationally 
demanding, and drives us to work with low resolution 
maps. Our goal is to determine the polarized CMB signal 
at large angular scale, so this does not limit the analysis. 

2.2.2. Sampling the index vector 

For fixed A, we sample from the conditional distribu- 
tion 

p(/3|A,d)ocp(d|/3,A)p(/3), (19) 

with prior probability p(/3). An analytic sample cannot 
be drawn from this distribution because the spectral in- 
dices are non-linear parameters. However, for a small 
number of parameters it is feasible to draw samples us- 
ing the Metropolis-Hastings algorithm. This algorithm 
has been described extensively in t he cosmolog i cal pa - 
rameter e stimation literature (e.g., [Knox et al.1 (|2QQll ): 
iLewis fc Bridle (2002); Dunkle^LeVaD (|2005f )). The sam- 
pling goes as follows. For each index parameter in turn, a 
trial step /3t is drawn using a Gaussian proposal distribu- 
tion of width ctt centered on the current /3 vector. Next, 
the current and trial /3 vectors are used to construct 
model vectors at each frequency, m^y = ^7 ock^i/ Az.. The 
current and trial posterior are then computed using 

- 2 lnp(/3| A, d) = - 2 lnp(/3), (20) 

with given in Eqn [5l The ratio of the trial to cur- 
rent posterior, r, is used to determine whether to move 
to the trial position (with probability r), or to stay at 
the original position (with probability 1 — r). This use 
of the Metropolis algorithm to draw a subset of parame- 
ters is commonly k nown as Metropolis-within-Gibbs (e.g. 
iGeweke fc Taniza ki (2001)), and has been u sed in a stron- 
omy to estimate Cepheid distances (Barnes et al. 2003). 
Other approaches to samplin g spectral index para meters 
have been considered in e.g.. lEriksen et al] (|2007l ). 

In regions of low signal-to-noise, it is necessary to im- 
pose a prior on the synchrotron spectral index, other- 
wise it is unconstrained and could take the 'fiat' in- 
dex of the CMB component, opening up large degen- 
eracies. We choose a Gaussian prior of —3.0 ± 0.3, mo- 
tivated by understanding of the synchrotron emission 
(|Rvbicki fc Lightmanlll979l ) and allowing for variations 
of the size ob s erved in the synchrotron intensity (e.g., 
iBennett et al.l (120031 )). This is combined with a uniform 
prior on the CMB and synchrotron amplitudes, and a 
Gaussian prior ± 0.21^ on the dust Q and U ampli- 
tudes. This parameterization and choice of prior does 
not guarantee that the marginalized means for the A 
and P param eters will be unbiase d estimators, as dis- 
cussed in e.g. lEriksen et al.l (|2007[ ). In the limit of low 
signal-to-noise, there is a larger volume of models with a 
fiatter (i.e. /3 tending to 0) synchrotron spectrum, allow- 
ing large CMB and synchrotron amplitudes of opposite 
sign. One approach is to modify the prior on the spectral 
indices to include an additional 'phase-space' factor. We 
discuss this further in ^ 

2.3. Processing the sampled distribution 

We form maps of each component from the mean of 
the marginalized distribution, 

. no 

(A,) = / p{Ak\d)AkdAk = — Va^„ (21) 



where the sum is over all tiq elements in the chain, and 
A^ is the zth chain element of the kih. component map. 
The covariance matrices for A/., including off-diagonal 
terms, are estimated using the same method, summing 
over the chain components. As an example, the covari- 
ance Cxy^k between pixels x and y for component k is 
computed using 

Cxy,k — {Ax^k^y,k) ~ {^x,k){^y,k) (22) 
-| riG 

= ^ T.«k - {A.,k)){Al^, - {Ay,k)). (23) 

For the synchrotron and dust components, we compute 
only the diagonal elements of the covariance matrices. 

2.3.1. Power in the E-mode 

To quantify the polarization anisotropy present in the 
Q and U maps, we use the coordinate-independent scalar 
and pseudo-scalar E and B modes commonly used in 
cosmological analysis (jSeljakI 119971 : iKamionkowski et al.l 
119971 ). Both polarization modes probe the evolution of 
the decoupling and reionization epochs and are gener- 
ated by Thomson scattering of a quadrupolar radiation 
pattern by free electrons. The anisotropy is quantified 
using the Cj^ ^ Cf^^ power spectra, where 

= {aLal:). (24) 

The spin-2 decomposition of the polarization maps, ctf^^ 
is related to the Q and U maps by 

[Q±iU]{x)=Y^ Yl T2(^irnT2y£m{x) (25) 

where ±2Ciim = cif^ ± iaf^ (Zaldarriaga fc Sell aki 1 19971 ). 

The first stars reionize the universe at redshift z^, pro- 
ducing a signal in the E-mode power spectrum propor- 
tional to r^, where r is the optical depth to reioniza- 
tion. In this analysis we use the approx imation to the 
optical depth used in iPage et al] (|2007lf ). estimated by 
varying only r and the power spectrum amplitude, such 
that the temperature anisotropy power at £ = 220 is 
held constant. Other cosmological parameters are fixed 
to fiducial values and the exact likelihood is computed 
as a function of r. The hkeliho od, ^(C^ld), given in Ap- 
pendix D of iPage et al] (|2007l ) , is evaluated using the 
marginal posterior mean Q/U maps and their covari ance 
matrix , which are processed as described in Page et al.l 
(|2007f ) to account for the P06 Galactic mask. 

2.4. Testing convergence 

The convergence of the chain is determi ned by applying 
the sp ectral convergence test described in lDunklev et al.l 
(|2005f ) to the spectral index parameters and a random 
subset of the amplitude parameters. This tests for con- 
vergence of the mean, but convergence of the covari- 
ance matrix is also required for estimating the power 
in the mapx at large scales. We check this by applying 
a jackknife test to the derived optical depth parameter 
r. For a chain that samples only the amplitude vec- 
tors, and has diagonal noise, 1000 iterations are typi- 
cally sufficient. With the Metropolis sampling step in- 
cluded, about 10,000 iterations are required, and when 
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off-diagonal noise is included, typically 50,000 iterations 
are necessary. 

3. SIMULATED POLARIZATION MAPS 

We simulate all-sky signal and noise maps at HEALPix 
^side=16, with Np = 3072 pixels, for the Stokes Q and 
U parameters, at the five WMAP frequencies (22.8, 33, 
40.7, 60.8, and 93.5 GHz). In the notation of ^ the data 
are modeled as 

= /(r.)Ai + {ly/iyKf' A2 + {yjvwf^K^, (26) 

where is a vector of length 2A/p containing the to- 
tal Q and U signal in antenna temperature at frequency 
V. The CMB signal map, Ai, is generated for a fidu- 
cial cosmological model by drawing multipoles aim from 
the theoretical power spectrum Q, and transforming to 
map space, using the 'synfast' routine in HEALPix. We 
choose a model with r = 0.1, with other cosmological pa- 
rameters give n by the best- fit three-year WMAP ACDM 
model (Spergel et al. '2007). For the synchrotron emis- 
sion, with amplitude map A2 defined at 23 GHz, we use 
the three- yea r WMAP Q and U low-resolution K-band 
maps (Page et" al.l 12007) . The frequency dependence is 
assumed to be a power law with = —3.0 over the full 
sky. For the thermal dust map, A3, defined at W band 
(94 GHz), we construct maps that are 5% polarized, us- 
ing the dust intensity template from IRAS, ex trapo- 
lated to 94 GHz using Model 8 of lFinkbeiner et all (l999 ) 
and degraded to A^side=16. We form 

Q3 = 0.051rfcos(27) (27) 
U3 = 0.05Irfsin(27), (28) 

where the polarization angles 7 are computed from 

the starli g ht po larization template described in 

iPage et all (|2007|). usi n g ob servations from iHeilesI 
(j2000r ): iBerdvudn et all (|200lL 120041 ). We assume a 
power law index with = 1.7 over the whole sky. This 
is a simplified model of the true sky, as we would expect 
the emission processes to result in some spatial and 
frequency dependence of the spectral indices. A realistic 
Galactic field model would also lead to a spatially 
dependent suppression of the dust polarization fraction 
(jPage et al.ll2007l : iMiville-Deschenes et al.ll2008[ ). 

The simulated maps are formed at each frequency v us- 
ing d^y = s^y + njy, where is a realization of the WMAP 
noise. We consider both diagonal noise realizations, and 
full noise realizations including pixel-pixel and Q-U cor- 
relations, drawn from the WMAP A^side=16 noise ma- 
trices. These corres pond to maps co-added by year and 
DA, as described in'Jarosik et al.l ((2007i) : iHinshaw et al.l 
(f2008) . The low resolution inverse noise matrices are not 
defined inside the processing mask, so we set the noise 
in these pixels to be large. Realizations are computed at 
each channel using SVD decompositions of the WMAP 
inverse noise matrices N~^, in antenna temperature. We 
then create A^side=8 simulations, with 768 pixels, by de- 
grading the A^side=16 simulated maps and inverse 
noise matrices using inverse weighting. 

We perform the Gibbs sampling using the four WMAP 
frequency bands K, Ka, Q, V (hereafter KKaQV). We 
do not include W band as stand ard, due to concern s 
about potential systematic effects ('Hinshaw et al.""2008). 
For the diagonal noise case this requires of order 10,000 
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Fig. 1. — Left: For maps with A^p = 768 pixels and A^^ (marked 
N in figure) synchrotron spectral indices with Gaussian priors of 
3.0±0.3 (dashed curve), the recovered index estimates are 
larger than the prior (solid curves) for simulations with no signal. 
The estimated index mean increases with the number of pixels 
sharing the same index. Right: The spectral index priors that are 
imposed (colored curves), in addition to ^5 = —3.0 ±0.3, to account 
for the volume effect. 

iterations for convergence. The chains are processed as 
described in §2.3l to obtain marginal posterior mean maps 
and errors for Ai (containing Q and U for CMB), A2 
(synchrotron), A3 (thermal dust) and /3 (synchrotron 
spectral index). 

3.1. Spectral index prior 

For the simulated maps with /3 = —3.0 and a Gaussian 
prior on the indices of —3.0 ± 0.3, we find the recovered 
marginal posteriors p{P\d) to be biased estimators. To 
explore this further, we sample the joint distribution for 
a noise-only simulation, drawn from the WMAP five-year 
noise maps. The resulting marginalized distributions for 
the Ni = 48 spectral indices are close to Gaussian, but 
centered on /3 = —2.6. The Gaussian distribution best 
fitting the samples is shown in the left panel of Figure 
[TJ When the number of spectral index pixels is increased 
to 96, and 192, the recovered distribution tends toward 
the input value of /3 = —3.0. This effect arises for our 
Galactic emission model, as there is a larger volume of 
phase space for shallow indices far from the prior central 
value. This is not a significant effect for individual pixels, 
but the phase-space volume relative to that of the prior 
central value increases for more pixels sharing a common 
index parameter. 

In the case of a single pixe l, a conimon approach 
is to adopt the Jeffreys' prior (|Jeffrevsl[T96lh . a non- 
informative prior distribution that is proportional to the 
square root of t he Fisher inform ation, as described for 
example in lEriksen et al.l (|2007l ). Imposing the Jeffreys' 
prior on the spectral index parameters in this model 
would give less weight to steep indices than shallow 
ones, accounting for their smaller effect on the likelihood. 
However, this would not overcome the phase space effect 
that arises from the increase in volume of models with 
shallow indices. 

While there may be alternative parameterizations that 
avoid this problem, we choose instead to modify the spec- 
tral index prior. We use p{/3) = Ps{/3)/ Lo{/3)^ where 
Ps is the Gaussian prior —3.0 ± 0.3, and I/Lq is an 
additional 'phase-space' factor. We estimate this fac- 
tor by evaluating the marginalized posterior, p{P\d), for 
a noise-only simulation with prior Ps{P)^ and defining 
Lo{/3) = p{/3\d)/ps{/3). The assumption is that Lq is 
an approximate description of the available phase-space 
volume. The factors, I/Lq, are modeled as Gaussian dis- 
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Fig. 2. — Comparison of input (left) and output (middle) component maps for a simulation with A^side 
P = ^Q2 ^ f/2 ig shown for the CMB (top), synchrotron at 23 GHz (middle), and dust at 94 GHz (bottom) 

the difference in standard deviations per pixel for the Q Stokes parameter. The dust component has a Gaussian prior on the dust Stokes 
parameters of [Q, U](n) = dz 0.2/d(n), where /^(n) is the FDS dust intensity (see text), which reduces the deviation per pixel and leads 
to the structure in the difference map. 



tributions and are shown in the right panel of Figure 
□ for A^i=48, 92 and 192, with (3 = -3.7 ± 0.5 for the 
Ni = 48 adopted in this analysis. These distributions 
are not exact, particularly at P < —4, but we check that 
the approximation is sufficient by resampling the distri- 
butions with the new prior, finding the recovered mean 
to he /3 = —3.0. We repeat the test for a Gaussian prior 
of —2.7 ±0.3, and further test that this prior only weakly 
depends on the inclusion of W band. 

This is not the optimal solution, as the prior is informa- 
tive in the presence of a signal and could therefore lead to 
parameter bias, in particular in the weak signal regime. 
We test simulations for two different signal maps with 
synchrotron index —3.0 and —2.7, using Gaussian priors 
Ps of —3.0 ±0.3 and —2.7 ±0.3 combined with the phase- 
space factor. As expected, when the prior ps matches the 
signal, the estimated index values are all consistent with 
the input signal. When it does not match the signal, an 
unbiased marginal index value is estimated in the highest 
signal-to-noise pixels, and the low signal-to- noise pixels 
tend to the central value of Ps. This is the behavior we 
expect, but does not confirm that our final estimated pa- 
rameters will be unbiased. However, since we have good 
astrophysical motivation to require that the synchrotron 
index is close to /3 = —3, we do not expect this to be a 
significant effect. We address this issue in this analysis 
by checking that our main conclusions are insensitive to 
the detailed choice of the spectral index prior, and defer 
an investigation of a more optimal parameterization to a 
future analysis. 



3.2. Simulations with diagonal noise 

We now apply the method to simulations with signal 
and diagonal noise properties. Figure [2] shows the input 
and output polarization amplitude maps for the CMB, 
synchrotron, and dust components, derived from the Q 
and U maps, P = \/Q^-U^, for the fiducial simulation 
with synchrotron index /3 = —3.0. The third column 
shows the difference in standard deviations per pixel for 
the Q Stokes parameter, S = (Qin — Qout)/o'Q. For the 
CMB maps the absolute difference between input and 
output maps is greatest in the Galactic plane where the 
foreground signal is high. However, since these effects are 
captured in the errors in the map, the deviation maps do 
not have a strong spatial dependence. The foi" the 
map, {A^ - A?^*)^C^^(A\" - A?^*), is 1270 for 1536 
pixels. This gives x^/pi^el < 1, due to the prior on the 
dust amplitude. When the prior on the dust amplitude is 
removed, the goodness of fit of the recovered CMB maps 
is — 1497 (x^ /pixel = 0.98), but the errors are sig- 
nificantly infiated. The synchrotron maps are recovered 
with /pixel = 1.04 for the Q and U maps. The dust 
maps are recovered with x^/pixel = 0.24, with struc- 
ture apparent in the deviation map. This is due to the 
prior: far from the plane the signal-to- noise is low and 
the dust tends to the prior central value of zero. Re- 
moving the prior on the dust amplitude, the goodness of 
fit of the recovered dust maps is x^/pixel = 0.97, close 
to 1 as expected. The total x^ of the estimated model, 
X;^[d^-m^]^N-^[d^-m^], is 2.087Vp without the dust 
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Fig. 3. — The distribution of the optical depth to reionization 
T for simulations with five-year WMAP noise levels. The optical 
depth is recovered for an individual simulation with r = 0.1 (black 
dashed line), for KKaQV and KKaQVW, and a r = simulation is 
consistent with zero power. Changing the index priors to (3^ = 2, 
or = — 2.7±0.3 have negligible effects on the recovered CMB 
power. Incorrectly modeling the synchrotron as a power-law for a 
simulation with an index that increases by c = 0.5 between 23 and 
94 GHz increases r by ~lcr. 

prior, and 3.51Np with the dust prior. This is consistent 
with the number of degrees of freedom {2Np — Ni in the 
absence of priors). 

The optical depth computed from this simulation is 
r = 0.091 ± 0.019 outside the P06 Galactic mask, with 
distribution shown in Figure [3l We test for bias by gen- 
erating 10 further diagonal noise and signal realizations 
of the model with r = 0.1, and find ensemble average of 
r = 0.098, consistent with the input but limited by small 
number statistics. We also test a simulation for KKaQV 
with no polarized CMB component. The recovered opti- 
cal depth to be consistent with zero, shown in Figure [3l 
Adding W band (KKaQVW) we find r = 0.098 ± 0.017. 

In this ideal case the simulation matches the parame- 
terized model, but in a realistic scenario the model will 
not perfectly describe the sky. For a Gaussian prior on 
the spectral index of /^^ — 2.7 ± 0.3, we find a negligible 
effect on r, with r = 0.094 ± 0.020. A larger increase of 
~l<j in r, to 0.111 ± 0.019, is seen when the spectral in- 
dex of the simulation shallows from -3.0 to -2.5 between 
K and W band, but is modeled by a pure power law. In 
this case the model does not remove enough synchrotron 
at high frequencies, leaving excess power in the CMB. 
For the thermal dust component we fix = 1.7 for all 
pixels in the fiducial case. Changing this to = 2 has a 
negligible effect on our results. Removing the dust prior 
altogether opens up long degeneracies between the dust 
and the CMB amplitudes, highlighting the importance of 
the dust intensity map to limit the polarized dust contri- 
bution. Modest changes of the dust polarization fraction 
prior to 15% or 25% have little effect on the estimated 
CMB signal. 

3.3. Simulation with full noise matrix 

When the full off-diagonal noise matrix is included we 
find of order 50,000 iterations are necessary for conver- 
gence. The foi" the map computed using the full co- 
variance matrix, {A^ - A?^*)^C5"^(A\" - A?^*) = 1220 
for 1536 pixels, consistent with the diagonal noise case 
(x^ = 1270). The recovered spectral index distributions 
are consistent with the simulated values, and the esti- 
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Fig. 4. — Marginalized 68% and 95% confidence levels for a sub- 
set of the A and (3 parameters for simulated maps (inputs shown 
dashed). The top and middle panels show the correlation between 
the CMB, synchrotron and dust Q and U amplitudes (in /iX) for an 
arbitrary single pixel (pixel 356 out of 768 using HEALPix nested 
ordering). The bottom panels show the correlation between Q and 
U for a single pixel (p439), and between two adjacent pixels (plO 
and pll). By marginalizing, rather than finding the maximum 
likelihood, the error on the CMB amplitude is inflated to account 
for foreground uncertainty. 

mated optical depth of a single simulation with input 
r = 0.1 is 0.110± 0.020. Testing a large set of simula- 
tions with full inverse noise matrices is beyond the scope 
of this analysis. To demonstrate the effect on the esti- 
mated CMB maps of marginalizing, Figure [4] shows two- 
dimensional marginalized distributions for a subset of pa- 
rameters, for a single pixel and between pixels. The top 
and middle rows show the correlation between the CMB, 
synchrotron, and dust Q and U components within a 
single pixel. The one-dimensional error on the marginal- 
ized CMB Q and U amplitudes is larger than the error 
obtained if the foreground amplitudes are fixed at their 
maximum likelihood amplitudes. The bottom row shows 
correlations between Q and U components within a pixel, 
and between adjacent pixels. If the inter-pixel correla- 
tions are ignored, the maps are recovered with incorrect 
noise properties. 

4. RESULTS FROM WMAP DATA 

We apply the sampling method to the five-year WMAP 
data, using low resolution coadded A^side=16 maps and 
inverse noise matrices as the inputs, and then degrading 
to A/'side=8 as for the simulations. Figures [5] and [6] show 
maps of the mean values and la errors for the marginal- 
ized CMB, synchrotron, and dust amplitudes. 

4.1. CMB polarization 

We first consider the CMB results. The noise patterns 
for both Q and U in Figure [6] are consistent with what we 
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Fig. 5. — Low resolution polarized Q (top) and U (bottom) maps of the CMB, synchrotron, and dust emission, estimated from the 
five-year K, Ka, Q, and V band maps using Gibbs sampling. Pixels inside the processing mask are grey. The CMB maps (left panels, 
thermodynamic temperature) do not show significant Galactic contamination in the plane. The synchrotron amplitudes (center, antenna 
temperature), are defined at K-band (22.8 GHz), and are consistent with the total K-band maps, with high Q and U emission from the 
North Polar Spur, and high Q emission in the Galactic plane at longitude 110 ^ / ^ 170. The dust amplitudes (right, antenna temperature) 
are defined at W-band (93.5 GHz), and have a Gaussian prior on Q and U of ± 0.2/^ where Id is the dust intensity. 
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Fig. 6. — Estimated la errors on the low resolution maps of the CMB (left), synchrotron (center), and dust (right) Q and U components, as 
shown in Figure [5] The CMB errors are more fully described by a covariance matrix, including pixel-pixel correlation and Q/U correlation, 
so the maps can be used for cosmological analysis. The errors on the dust maps (right) are dominated by the prior that limits the dust 
polarization fraction to 20%. The middle panels clearly show the two sources of uncertainty in our CMB polarization maps: detector noise 
in the ecliptic plane (which traces a sideways S in the map) and foreground removal uncertainties in the galactic plane. 



expect: in regions of low Galactic emission, the errors are 
dominated by instrumental noise. As the Galactic plane 
is approached, errors from foreground uncertainty begin 
to dominate, in particular where the dust contribution is 
most uncertain. This is a real advantage of the method: 
rather than imposing masks, the method inflates the er- 
rors where the foregrounds are brightest. As opposed to 
template cleaning, which produces CMB maps at each 
frequency observed, this method recovers a single Q and 
U polarization map, and so has higher signal-to-noise 
than any of the individual template-cleaned maps. There 
is some indication of structure in the CMB signal, par- 
ticularly in the Q Stokes map at the Galactic anti-center, 
and in the U map in the region of the North Polar Spur. 
This is consistent with noise. Outside the P06 mask the 
maps are morphologically similar to the template-cleaned 
CMB maps co- added over Ka, Q, and V bands. The cor- 



relation coefficients for the pixels within the CMB maps 
have a maximum value of ~ 40%, with rms correlation 
of - 2%. 

We compare the power at each multipole, outside the 
P06 mask, to the template-cleaned case from the main 
WMAP analysis (Gold et al. 2008) in Figure^ Using the 
method described in INolta et al.l ()2008l ). at each multi- 
pole the conditional likelihoods is computed as a func- 
tion of Cf^^ with all other multipoles held fixed, for 
2 < i < 7. The results are consistent, although this 
analysis finds more power at ^ = 4 and 5. Computing 
the likelihood as a function of r we find r = 0.090 ±0.019 
for the Gibbs-sampled maps outside the P06 mask, which 
is consistent with the results obtained through template 
cleaning, which give r = 0.086 ± 0.016 for the KaQV 
data combination. Obtained using a different methodol- 
ogy and accounting for foreground marginalization, this 
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Fig. 7. — Conditional likelihoods for the CMB EE multipole mo- 
ments estimated from the polarization maps described in this anal- 
ysis (black) , compared to the template cleaned maps described in 
IGold et al.l (pOOl ) (red). They are computed as in Nolta et alj 
(|2008l ) by fixing all other values at the fiducial ACDM values. 
There is agreement between the two methods. 

adds confidence in the detection of the CMB E-mode 
polarization signal. The spatially varying synchrotron 
index appears not to cause a significant difference, as we 
obtain a similar mean when the synchrotron spectral co- 
efficients are fixed at the best-fit values found in the tem- 
plate cleaning. This also confirms that the choice of the 
spectral index prior is not affecting the estimated CMB 
power in the maps. We find similar limits on the opti- 
cal depth when W band is included (KKaQVW), with 
r = 0.085 ± 0.017. However, we choose not to use W 
band in the standard analysis, as discussed earlier. The 
Gaussian prior on the dust Q and U parameters does af- 
fect the CMB signal constraint: removing it significantly 
weakens the limit on the large-scale power, but changing 
it to e.g., 15% or 25% has little effect, consistent with 
tests on simulations. A significantly tighter prior is not 
physically motivated, and could lead to bias in the re- 
covered signal. 

4.2. Synchrotron polarization 

The synchrotron maps shown in Figure [5] are similar 
to the total K-band maps (Hinshaw et al. 2008). The 
difference between the estimated synchrotron amplitude, 
and the K-band amplitude, is < 5 fiK outside the P06 
mask, and < 8 jaK in the Gal actic plane. As ob served in 
the three-year WMAP data (iPage et al.ll2007f ). the sig- 
nal is dominated by emission from the No rth Polar Spur, 
marked on the microwave sky map in Hins haw et al I 
([2007), as well as what is often known as the 'Fan re- 
gion' (e.g.. IWolleben et al.l (|2006l )). centered on Galactic 
coordinates / ~ 140, 6^5. The synchrotron emission 
dominates the signal at low frequencies, and so the un- 
certainty in the synchrotron Q and U maps, shown in 
Figure [H is dominated by instrument noise, with only 
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Fig. 8. — The estimated synchrotron spectral index (top), and 
la errors (bottom), estimated in 48 HEALPix A^gi(^g=2 pixels. A 
Gaussian prior of (3s = —3.0 zb 0.3 is imposed in each pixel. In 
regions of low signal-to-noise (near the ecliptic poles), the prior 
drives the spectral index estimate, so we mask the index for pixels 
with cr(/3) > 0.25. The mean index averaged over the unmasked 
pixels is -3.03 ± 0.04 for prior -3.0 ± 0.3. 



a small contribution from marginalization in the Galac- 
tic plane. Figure [8] shows the mean synchrotron spectral 
index estimated in the 48 pixels, together with la er- 
rors. In regions of low synchrotron signal-to-noise the 
index is driven by the prior, so we mask pixels with 
cr(/3) > 0.25. There is a preference in the North Po- 
lar Spur and Fan region for an index of ~ —3: the index 
averaged over regions with a{P) < 0.25 is —3.03 ± 0.04. 
The estimated indices of these pixels are also shown in 
Figure [9l ordered by increasing errors to highlight the 
high signal-to- noise pixels on the left. Cutting the sky 
into high and low latitude at 6 = 20 degrees, the low 
latitude weighted mean, —3.00 ± 0.04, is consistent with 
the high latitude —3.08 ± 0.06. This contrasts with the 
more signif icant steepening o f the index with latitude ob- 
served in Kogut et al.l (|2007[ ) in the analysis of the three- 
year WMAP polarization data, although here we exclude 
part of the plane. These results coming from alternative 
analyses will be more easily assessed with higher signal- 
to-noise data. We check that the latitude dependence is 
unaffected when the phase-space factor in the spectral 
index prior, a Gaussian — 3.7 ± 0.5, is modified to have 
width 1.0 or mean —4.0. 

An extrapolation of the Haslam 408 MHz synchrotron 
intensity maps to WMAP frequencies suggests that a 
shallower index is preferred in the plane, compared to the 
polarization, in the absence of anomalous dust. However, 
our estimated polarization index is consistent with the 
synchrotron intensity maps derived by the MEM method 
()Gold et al.l 1200 8). Degrading the synchrotron compo- 
nent map in each band to A^side=2, and fitting a power 
law to the data in K, Ka, Q, and V band for each pixel, 
the spectral index values are best-fit by a power-law with 
index /3 = —3.1 for the pixels at 6 = 0. 

We check the estimate of the spectral indices by re- 
running the analysis with a Gaussian synchrotron index 
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Fig. 9. — Marginalized mean and la error bars for the syn- 
chrotron spectral index parameters with highest signal-to-noise for 
the five-year WMAP data, for Gaussian priors of — 3.0 ± 0.3 and 
—2.7 =b 0.3. The pixels are ordered with increasing Icr errors, with 
the highest signal-to-noise pixels on the left. The prior has only a 
small effect on the estimated indices. 

prior, Ps{/3)^ of —2.7 ± 0.3 in addition to the phase-space 
factor. The estimated indices are compared pixel by pixel 
in Figure [9l The prior has only a small effect on the 
high signal-to-noise index estimates. The (j{/3) < 0.25 
pixel average in this case is —2.94 ±0.04. The estimated 
CMB power is little affected by this change in prior, with 
r = 0.092 ±0.019. 

4.3. Dust polarization 

The dust polarization map has a low signal-to-noise ra- 
tio, particularly far from the plane, as we only fit data in 
the K-V bands. In these regions the prior dominates the 
estimate of the dust amplitude, making it hard to draw 
conclusions about the dust component. The error in Q 
and U is driven by the prior on the polarization ampli- 
tudes and so is morphologically identical to the FDS dust 
intensity map. This explains why the error far from the 
plane is low even though the dust is only poorly mea- 
sured. The fractional polarization outside P06 is typ- 
ically only 1-2%, where we use the degraded FDS dust 
map to trace intensity, and rises to ~10% in some regions 
of the plane. This is lower than the ~ 4% estimated 



in (jKogut et all 120071 : IGold et al.ll2008[ ). However, these 
maps are only estimated for z/ < 61 GHz, and in regions 
of low dust the prior prefers zero polarization. The frac- 
tional polarization estimate also assumes that the FDS 
map accurately traces the dust intensity in the WMAP 
frequency range. Inclusion of higher frequency data will 
allow us to learn more about the polarized dust. 

5. CONCLUSION 

We have used Bayesian parameter estimation to es- 
timate low resolution polarized CMB maps, marginal- 
ized over foreground contamination. These may then 
be used as inputs for a likelihood analysis. The emis- 
sion model is parameterized accounting for physical un- 
derstanding of the Galactic emission. The method has 
been tested on simulated maps, and found to produce 
unbiased estimates for the CMB power, quantified by 
the optical depth to reionization. With the five-year 
WMAP data we find a consistent result compared to tem- 
plate cleaning, with r = 0.090 ± 0.019 from this method 
and 0.086 ± 0.016 from the standard template-cleaning 
method. This method captures the increase in errors 
where foreground uncertainty is larger, so depends less 
on a Galactic mask. Estimates of the polarized Galac- 
tic components indicate a synchrotron spectral index of 
order /3 = —3.0 in the Fan region in the Galactic anti- 
center, and the North Polar Spur area. 
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